********************************************************************************
* make a skinny version of the data
********************************************************************************

use dataset_voters, clear
local X anywells highfracking zipcode year turnout anid
keep `X'
save dataset_voters_tmp, replace

********************************************************************************
* set up data set to store results
********************************************************************************

clear
set obs 100
gen B = .
gen Bhf = .

forvalues k = 1(1)100 {

di `k'
di  "`c(current_time)'"
qui{

preserve

********************************************************************************
* make the placebo fracking variables
********************************************************************************

use dataset_voters_tmp, clear
collapse anywells highfracking, by(zipcode year)
gen placebo = 0
gen placebohf = 0 if highf == 1
encode zipcode, gen(id)
xtset id year

* make it for all zip codes

forvalues i = 2004(4)2016 {
	su anywells if year == `i'
	local b`i' = r(mean)
}	
	
replace placebo = rbinomial(1, `b2004') if year == 2004
replace placebo = l4.placebo if year == 2008
replace placebo = l4.placebo if year == 2012
replace placebo = l4.placebo if year == 2016

local myb = `b2008' - `b2004'
replace placebo = placebo + rbinomial(1, `myb') if placebo == 0 & year == 2008
replace placebo = l4.placebo if year == 2012
replace placebo = l4.placebo if year == 2016

local myb = `b2012' - `b2008'
replace placebo = placebo + rbinomial(1, `myb') if placebo == 0 & year == 2012
replace placebo = l4.placebo if year == 2016

local myb = `b2016' - `b2012'
replace placebo = placebo + rbinomial(1, `myb') if placebo == 0 & year == 2016

* make it for high-fracking zip codes

forvalues i = 2004(4)2016 {
	su anywells if year == `i' & highfracking == 1
	local b`i' = r(mean)
}

replace placebohf = rbinomial(1, `b2004') if year == 2004 & highf == 1
replace placebohf = l4.placebohf if year == 2008 & highf == 1
replace placebohf = l4.placebohf if year == 2012 & highf == 1
replace placebohf = l4.placebohf if year == 2016 & highf == 1

local myb = `b2008' - `b2004'
replace placebohf = placebohf + rbinomial(1, `myb') if ///
	placebohf == 0 & year == 2008 & highf == 1
replace placebohf = l4.placebohf if year == 2012 & highf == 1
replace placebohf = l4.placebohf if year == 2016 & highf == 1

local myb = `b2012' - `b2008'
replace placebohf = placebohf + rbinomial(1, `myb') if ///
	placebohf == 0 & year == 2012 & highf == 1
replace placebohf = l4.placebohf if year == 2016 & highf == 1

local myb = `b2016' - `b2012'
replace placebohf = placebohf + rbinomial(1, `myb') if ///
	placebohf == 0 & year == 2016 & highf == 1

********************************************************************************
* merge back to individual data and do the regressions
********************************************************************************

merge 1:m zipcode year using dataset_voters_tmp
reg turnout placebo i.year, a(anid)
local B = _b[placebo]
reg turnout placebohf i.year, a(anid) 
local Bhf = _b[placebohf]
 
restore

replace B = `B' if _n == `k'
replace Bhf = `Bhf' if _n == `k'

* save

outsheet using "results.csv", comma replace

}

}

********************************************************************************
* graph
********************************************************************************

* get observed point estimates
/*
preserve
use dataset_voters, clear
areg turnout anywells i.year, a(anid) cl(zipcode)
areg turnout anywells i.year if highfracking == 1, a(anid)
*/

preserve

insheet using "results.csv", clear 

count if abs(b) >= 0.0154

#delimit;

gr tw
	(hist b, fcol(white) freq)
	,
		xline(-0.0154 0.0154, lpat(dash) lwid(thick)) 
		xlab(-0.02(0.01)0.02)
		ylab(, angle(horiz))
		plotregion(style(none))
		xtitle("Placebo estimate")
		title("All states")
		subtitle(`r(N)' out of 100)
		text(20 -0.008 "{&larr} Observed estimate")
		name(g1, replace)
		;

#delimit cr

count if abs(bhf) >= 0.0195

#delimit;

gr tw
	(hist bhf, fcol(white) freq)
	,
		xline(-0.0195 0.0195, lpat(dash) lwid(thick)) 
		xlab(-0.02(0.01)0.02)
		ylab(, angle(horiz))
		plotregion(style(none))
		xtitle("Placebo estimate")
		title("High-fracking states")
		subtitle(`r(N)' out of 100)
		name(g2, replace)
		;
		
#delimit cr

gr combine g1 g2, xsize(2) ysize(1)

gr export "_output/figureA5.pdf", replace

erase dataset_voters_tmp.dta

********************************************************************************
* end
********************************************************************************
